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Abstract: We present a wave function approach to study the evolution of a 
small system when it is coupled to a large reservoir. Fluctuations and dissipation 
originate in this approach from quantum jumps occurring randomly during the 
time evolution of the system. This approach can be applied to a wide class of 
relaxation operators in the Markovian regime, and it is equivalent to the stan- 
dard Master Equation approach. 

Published in AIP Conference Proceedings 275, Thirteenth International Con- 
ference on Atomic Physics, Munich, Germany 1992; Editors: H. Walther, T. W. 
Hansch, and B. Neizert. 



The problem of dissipation plays a central role in Atomic Physics and Quan- 
tum Optics. The simplest example is the phenomenon of spontaneous emission, 
where the coupling between an atom and the ensemble of modes of the quantized 
electromagnetic field gives a finite lifetime to all excited atomic levels. Usually 
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the dissipative coupling between a small system and a large reservoir can be 
treated by a master equation approach [H [21 El H]; one writes a linear equa- 
tion for the time evolution of the reduced system density matrix p s = Tr rea (p), 
trace over the reservoir variables of the total density matrix. If we denote the 
hamiltonian for the system H$, this equation can be written: 

PS = t\PS, H S ] + £ rc lax(Ps) • (1) 

In (CQ), £ re iax is the relaxation superoperator, acting on the density operator p$. 
It is assumed here to be local in time, which means that ps(t) depends only 
on ps at the same time (Markov approximation). All the system dynamics can 
be deduced from (JT|). One can calculate one time average values of a system 
operator A: a(t) = (A)(t) = Ti(ps{t)A), and also, using the quantum regression 
theorem [5], multi-time correlation functions such as (A(t + r)B(t)). 

We present here an alternative treatment based on a Monte-Carlo evolution of 
wave functions of the small system (MCWF) [6j [7J [H [9] . This evolution consists of 
two elements: evolution with a non hermitian hamiltonian, and randomly decided 
"quantum jumps", followed by wave function normalization. This approach, 
which is equivalent to the master equation treatment, has two main interests. 
First, if the relevant Hilbert space of the quantum system has a dimension N large 
compared to 1, the number of variables involved in a wave function treatment (~ 
TV) is much smaller than the one required for calculations with density matrices 
(~ iV 2 ). Second, new physical insight may be gained, in particular in the studies 
of the behavior of a single quantum system. 



1 The MCWF procedure 

The class of relaxation operators that we consider here is the following: 

Wps) = -\ £ (Clomps + PsClC m ) + J2 C mPs Cl . (2) 

This type of relaxation operators is very general and it is found in most of the 
Quantum Optics problems involving dissipation. In (T5]), the C TO 's are operators 
acting in the space of the small system. Depending on the nature of the problem 
there can be one, a few or an infinity of these operators. 

For the particular case of spontaneous emission by a two-level system with 
one stable ground state g and one excited state e with a lifetime T -1 , there is 
just a single operator C\ = VT\g)(e\ in the relaxation operator ([2]), and one can 
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check that (j2J) indeed leads to the well known relaxation part of the optical Bloch 
equations: 

f (Ps) ee = -T( PS )ee { (Ps) eg = -(T/2)( PS ) eg 

\(PS)„ = ^ee \(ps) ge = -(T/2)( PS ) ge . 

We now present the procedure for evolving wave functions of the small system. 
Consider at time t that the system is in a state with the normalized wave function 
\4>(t)). In order to get the wave function at time t + 5t, we proceed in two steps: 

1. First we calculate the wave function + St)) obtained by evolving 

\4>{t)) with the non hermitian Hamiltonian: 

H = H s -^Y,CLC m . (4) 

This gives for sufficiently small 5t: 

mt+st)) = (i-^y<f>(t))- (5) 

Since H is not hermitian, this new wave function is clearly not normalized. 
The square of its norm is: 



(^\t + 8t)\^(t + 5t)) = <<K*)l(l + ^)(l-^ 



= l-8p (6) 

where Sp reads: 

5p = 5t l -{<p{t)\H-H^{t)) = H 5 Pm (7) 

5p m = 5t (<p(t)\ClC m \<P(t)) > 0. (8) 

The magnitude of the step 5t is adjusted so that this calculation at first 
order is valid; in particular it requires 8p <C 1. 

For the particular case of the two-level atom problem, the non hermitian 
Hamiltonian is 

H = H 8 ~\e)(e\. (9) 

This amounts to adding the imaginary term —ihT/2 to the energy of the 
unstable excited state, as usual in scattering theory. 
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Figure 1: The possible quantum jumps in the Monte-Carlo evolution 



2. The second step of the evolution of \<j>) between t and t + St consists in a 
possible "quantum jump" (Fig. 1). The various possible "directions" for 
the jumps are given by the C m operators, and the probability for making 
a jump in the "direction" of a particular C m is given by Sp m given in (jSJ). 
The new normalized wave function after such a jump is given by: 



with a probability Sp n 



m+st)) 



C m \<P{t)) 

\c m W))\\ 



(10) 



Using (Jj[J), we find that the total probability for making a jump is Sp. In 
the no-jump case, which occurs then with a probability 1 — Sp, we take as 
new normalized wave function at time t + St: 



with a probability 1 — Sp = 1 — V] Sp Tt 



m+st)) 



\^\t + St)) 

IW + ^)>ll ' 

(ii) 



Consider again as an example the particular case of the spontaneous emission 
of a two-level atom. The wave function at time t can be written as: 



m)) = a{t)\e) + (3{t)\9) 



(12) 
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Since there is a single C m operator in this case, there is only one possible type 
of quantum jump. The probability for this quantum jump is: 



and the wave function after the jump is simply \<f)(t+St)) = \g). If no jump occurs, 
the wave function at time t+St is similar to ([12]) . with the coefficients a(t+St) and 
(3(t + St) deduced from a(t) and /3(t) using the evolution with the non hermitian 
hamiltonian (Q. We see for this particular case that the Monte-Carlo evolution 
can be understood as the stochastic evolution of the atomic wave function if a 
continuous detection of the emitted photons is performed. The probability for 
detecting a photon during a particular time step St is indeed equal to Sp given in 
( TT3l) . and the new wave function after the detection, according to the standard 
quantum measurement theory, corresponds to the atom in its ground state g. 

It is actually quite a general result that the Monte-Carlo evolution outlined 
above represents a possible history of the system wave function with a suitable 
continuous detection process taking place [6J [8]. Although this procedure does 
not make any reference to measurements on the system, it may be useful, in 
order to get some physical understanding for the result of the simulation, to 
refer to such a continuous detection process, as if it was really performed. We 
note in this respect that one might possibly consider several different continu- 
ous detection processes for a given quantum system. The various sets of C m 's 
associated to each of these detection schemes can be deduced from each other 
by linear combinations, the relaxation equation (j2J) remaining then of course 
unchanged [?]. 

2 Equivalence with the Master Equation 

With this set of rules we can propagate a wave function | </>(£)) in time, and we 
now show that this procedure is equivalent to the master equation (1T1). More 
precisely we consider the quantity a(t) obtained by averaging a(t) = \<f>(t)) (<f>(t)\ 
over the various possible outcomes at time t of the MCWF evolutions all starting 
in |0(O)), and we prove that a(t) coincides with ps(t) at all times t, provided 
they coincide at t — 0. 

Consider a given MCWF \4>(t)) at time t. At time t + St, the average value 
of a(t + St) is: 



Sp = T\a\ 2 St 



(13) 



a(t + St) 
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h Vm \\c m \m)\\\\c m \m)\\ 1 ] 



which gives, using ([5j) , (IHJ) and ( TTOT) : 



i8t, 



a(t + 5t) = a(t) + -jr[(r(t), H s ] + St C rclax (a(t)). (15) 

We now average this equation over the possible values of a(t) and we obtain: 

-£ = l.[a,H s }+C IC Ua). (16) 

This equation is identical to the master equation ([1]). If we assume that ps(0) = 
|0(O)) (0(0) |, a(t) and ps(t) coincide at any time, which demonstrates the equiv- 
alence between the two points of view. In the case where ps{0) does not corre- 
spond to a pure state, one has first to decompose it as a statistical mixture of 
pure states, p(0) = J2Pi\Xi)(Xi\ an d then randomly choose the initial MCWFs 
among the with the probability law p^. 

As mentioned in the introduction, the master equation approach and the 
reduced density matrix give access to one time average values a(t) = (A)(t) = 
Tr(ps{t)A), which can now also be obtained with the MCWF method. One 
calculates, for several outcomes |0^(t)) of the MCWF evolution, the quantum 
average (0^(t)|v4|0W(t)), and one takes the mean value of this quantity over the 
various outcomes \(j)^\t)}: 

(A) (n) (t) = -j:(<P {l) (t)\A\^\t)) ■ (17) 

For n sufficiently large, f|T6|) implies that (A)r n ){t) — (A)(t). The ability to 
provide expectation values for any system operator makes the MCWF method a 
computational tool which may be much more efficient than the numerical solution 

of (l) [HE]. 

As an example of the agreement between the master equation approach and 
the MCWF approach, we have calculated by both methods the excited state 
population of a two-level atom coupled to a coherent laser field. The parameters 
for this Rabi nutation are a zero detuning 5 between the laser and atomic fre- 
quencies, and a Rabi frequency Q = 3T. In Fig. 2a, we show the excited state 
population for a single "history" for |0(t)). One finds a continuous evolution 
for this population oscillating between and 1, interrupted by random quantum 
jumps projecting the atomic wave function into the ground state. In Fig. 2b, we 
indicate the MCWF result obtained with the average of 100 wave functions. It 
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shows a damped oscillation as a result of the dephasing of the individual oscilla- 
tions due to the randomness of the various quantum jumps. The MCWF result 
is in good agreement with the one derived from the master equation (Optical 
Bloch Equations). Note that the purpose of this example is to illustrate the 
convergence of the two methods, and not to provide an efficient way of treating 
two-level atom problems. For such a small system, there is of course no gain in 
computing time by using the MCWF method instead of the master equation. 




time (unit: 1/r) 



Figure 2: (a) Time evolution of the excited-state population of a two-level atom 
in the MCWF approach. The dashed lines indicate the projection of the atomic 
wave function onto the ground state (quantum jump), (b) Excited state popula- 
tion averaged over 100 MCWF starting all in ground state at time 0. The dotted 
line represents the master equation result. 

It appears clearly that the equivalence of the Master Equation and MCWF 
approaches does not depend on the particular value of the time step 5t. From a 
practical point of view, the largest possible 5t is preferable, and one might benefit 
from using a generalization of (jHj) to a higher order in 5t, as for example a 4th 
order Runge-Kutta type calculation. The only requirement on 5t is that the 
various where the hrji are the eigenvalues of H, should be small compared 

to 1. Of course we assume here that those eigenvalues have been simplified as 
much as possible in order to eliminate the bare energies of the eigenstates of Hg- 
For instance, for a two-level atom with a transition frequency uja coupled to a 
laser field with frequency ljl, one makes the rotating wave approximation in the 
rotating frame so that the |^j|'s are of the order of the natural width T, the Rabi 
frequency Q or the detuning 5 = u>l — uja\ they are consequently much smaller 
than u>a ■ 

One might wonder whether there is a minimal size for the time step 5t. In 
the derivation presented above, it can be chosen arbitrarily small. However one 
should remember that the derivation of involves a coarse grain average of 
the real density operator evolution. The time step of this coarse grain average 
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has to be much larger than the correlation time r c of the reservoir, which is 
typically an optical period for the problem of spontaneous emission. Therefore 
one should be cautious when considering any result derived from this MCWF 
approach involving details with a time scale of the order of or shorter than r c , 
and only St larger than r c should be applied. This appears clearly if one starts 
directly from the interaction Hamiltonian between the system and the reservoir 
in order to generate the stochastic evolution for the system wave function [6]. 
The condition St ^> r c is then required to prevent quantum Zeno type effects 
|10j . This restriction is discussed in detail in [TTJ in connection with quantum 
measurement theory. 

3 MCWF and other stochastic approaches 

The problem of stochastic wave function evolution in connection with the treat- 
ment of dissipative systems in quantum optics has recently received a lot of 
attention. In the context of non classical field generation, Carmichael [E] has 
proposed an approach named "quantum trajectories", inspired by the theory 
of photoelectron counting sequences p2] and quite similar to the spirit of the 
present work. 

For simple atomic systems (2 or 3 levels) coupled to the electromagnetic field, 
the dynamics can be interpreted in terms of one or a few delay functions, which 
give the probability distribution of the time intervals between the emission of two 
successive photons [13j [14j [15]. When these functions are known analytically, 
they can generate a very efficient Monte-Carlo analysis of the process: just after 
the emission of the n th fluorescence photon at time t n , the atom is in its ground 
state and the choice of a single random number is sufficient to determine the 
time t n+ i of emission of the n + I s * photon. This type of Monte-Carlo analysis 
has been used in [16] to simulate an atomic beam cooling experiment, and in 
[H] to prove numerically the existence of dark periods in the fluorescence of a 3- 
level atom (quantum jumps). Very recently, laser cooling of atoms using velocity 
selective coherent population trapping [17J and lasing without inversion [18] have 
been analyzed by this type of Monte-Carlo method. 

Unfortunately, the delay function cannot be calculated analytically for com- 
plex systems involving a large number of levels. Nevertheless, it is possible 
to generate a Monte-Carlo solution for this problem in which a single random 
number determines the time of emission of each fluorescence photon [9]. The 
evolution of the system between two quantum jumps has to be integrated step 
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by step numerically, so that the amount of calculation involved is similar to the 
one required by the method presented in this paper. 

Stochastic approaches have also been introduced in the context of either stan- 
dard [Hi [20] or quantum non demolition [21] [22j [23] measurements of photon 
number in a given mode of the electromagnetic field. A sequence of random 
quantum jumps resulting from successive measurements asymptotically leads to 
a reduction of the field state into a Fock state \n), whose probability distribu- 
tion is equal to the initial photon number distribution for the case of the non 
demolition measurement. The main interest of these stochastic approaches, as 
compared with the usual master equation treatment, is to give explicit individual 
histories of the quantum field state in a measurement sequence. This is particu- 
larly valuable if one wants to optimize the measurement sequence in order to get 
a complete information on the field state with a minimum number of measure- 
ment processes |22j. On the other hand, these stochastic calculations still mostly 
deal with density matrices and their authors do not seem to consider them as 
more efficient ways of computing than the master equation. 

Another class of stochastic equations for system wave functions, which is also 
equivalent to the master equation ([1]), has been introduced by Gisin and Percival 
[21] (see also the work by Diosi |25j). In this approach, only continuous stochastic 
equations are considered. The complex Ito stochastic process is given by: 

\d<p) = - l -H s dt\<p) + Y,((Cl)C m -^ClC m -^(Cl)(C m ))\<P)dt 

+ E(c m -<c m ))w^f ( 18 ) 

where (C m ) = ((p\C m \4>), and where the d^ m are independent complex Wiener 
processes pi]: 

d£,m = 

WMmWQ = WUW&) = Sm,ndt (19) 

WdUWU = 

Carmichael has shown that for the particular case of the homodyne detection 
of the fluorescence light, the Quantum Jump formalism can be transformed into 
such a continuous stochastic equation [8j. Actually this proof can be extended 
to the most general case: the first step is to write the relaxation operator £ re iax 
as: 

CreUPs) = ~H {Dl, £ D m , £ p s + p s Dl >£ D m>£ ) + £ D m , £ p s Dl >£ (20) 

^ m,e m,e 
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where e — ±1 and where the D m £ are defined as: 

ull + eC m , . 

m ' £ = V2 (21) 

One easily shows that £ re iax in (1201) is identical with the one in (2). The coefficient 
fi is arbitrary at this stage; fi 2 has the dimension of the inverse of a time, and 
we just require in the following fi 2 ^> \rj\, where hrj is a typical eigenvalue for H 
(for the two-level atom case, r\ ~ r, Q,5). Using the set of operators -D m , e , we 
can now perform a Monte-Carlo evolution of the wave function, equivalent to the 
master equation (QQ). Because of the large magnitude of /i 2 , this simulation with 
the D m>£ operators involves much more quantum jumps in a given time interval 
At than a simulation done with the C m 's. But the change of the wave function 
in a given quantum jump: 

DmM (22) 



\\DrnMW 

is very small since D m e is nearly proportional to the identity operator t. In the 
limit of very large fi, the Monte- Carlo evolution of the wave function therefore 
tends towards a continuous stochastic evolution. In Carmichael's homodyne 
detection problem, the form (12T1) for the D m ^ £ has a clear interpretation. These 
jump operators correspond to the detection of a photon after one has mixed the 
light emitted by the atomic system with a local oscillator field. The parts in /ill 
and C m correspond respectively to the field originating from the local oscillator 
and the field emitted by the atom. The condition fi 2 3> \r)\ just states that 
the intensity of the local oscillator is much larger than the intensity of the light 
emitted by the atom, as usual in homodyne detection. 

To prove the equivalence with (fl8l) . we choose a time interval At such that 

fi- 2 < At < | ^7 1 1 . (23) 

This implies that the number of jumps N mj£ occurring with a given operator 
D m e during At will be large compared to 1 since fj 2 At 3> 1, but at the same 
time we expect only a small change in the system wave function since \r)\ At <C 1. 
The operator O describing the action of all those jumps during At is a product 



of the various D m:£ and it can be approximated at order 1 in y At \rj\ by: 

/ \ N , 

4= (u + -V(iV mi+ -iV m ,_) C m ) (24) 



where N = J2 m E ^m,e is the total number of jumps occurring during At. The 
wave function at time t + At can now be written before normalization: 

m + At)) = (l - % -H s M - \ At £ ClC m + £ Nm > + ~ Nm >- C m ) \m (25) 
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where we have taken into account both the non-hermit ian evolution during At 
and the effect of the multiple quantum jumps. The numbers of jumps N m ^ £ are 
poissonian random variables with an average value and a standard deviation 
given by: 



(26) 



AN m , £ ~ (27) 

where the average value (C m + C m ) is taken in \<fi(t)). In the limit of large N mi£ , 
we can approximate the random variable iV m + — iV m) _ appearing in (|25|) by: 

Nm ' + ~ ^= = At(C m + C*> + AC m (28) 

where A( m is a real gaussian random variable with zero mean and a standard 
deviation equal to y/At. Finally we normalize the wave function (!25j) and we 
obtain: 

|A0(t + At)> = - ~H s \<f>{t)) At 

+ \T. (<G» + cl)c m - clc m -\{c m + cl) 2 ) \4>{t))*t 

+ W^Cm-iQn + Cl^imMm. (29) 

m 

In (1291) . we have kept terms linear in A( m and At, and we have replaced all the 
quadratic terms A( m A( m > by their mean At 8 m , m i. In the limit fi — > +00, At — > 
0, this equation can be understood as a Ito stochastic equation, corresponding 
to a real version of (jl8j) . 

The exact form of ( TLB"]) can be recovered by taking a slightly more complicated 
set of D m e operators: 

D m , £ = 1X1 + 2 £ ° m with e = ±l,±i (30) 

and by performing an appropriate global phase change of the wave function. The 
continuous stochastic equation f|T8l) is therefore a limiting case of the quantum 
jump formalism presented here, and it also has an interpretation in terms of 
a detection scheme: the information concerning the system is mixed with a 
"classical field" /xl, and the sequence of quantum jumps deduced from the whole 
set of mixed components D m ^ £ allows one to determine the subsequent system 
evolution. Note that on the contrary, the jumps deduced from the action on the 
system wave function of a single mixed component, such as the -D mi+ 's, are not 
sufficient to determine this system evolution. 
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4 Conclusion 



We have presented a stochastic evolution for the wave function of a system cou- 
pled to a reservoir in the Markovian regime. Each time step in this stochastic 
evolution consists in two parts: a Hamiltonian but non hermitian evolution and 
a possible quantum jump. We have proved the equivalence of this Monte-Carlo 
Wave function approach with the master equation treatment. We have also 
shown that this simulation with Quantum jumps can be transformed into a con- 
tinuous stochastic evolution of the wave function, similar to the one of [24J . 

This approach provides a computational tool which is often more efficient 
than the standard master equation treatment for systems with a number of states 
N ^> 1 (for a detailed discussion see [7J). Indeed a wave function involves only 
iV components while a density matrix is described by iV 2 terms. This method 
has already been applied successfully to problems such as the study of the limits 
of laser cooling in 2 dimensions [26], or the calculation of the spectrum of the 
light emitted by an assembly of cold atoms [27] . Problems such as the study of 
collisions between cold atoms, or non linear mixing of quantum fields may also 
benefit from such an approach. 

We have emphasized that this simulation is in many practical cases directly 
connected to a measurement sequence performed on the system. Each Monte- 
Carlo trajectory is a possible history for the individual quantum system. In 
this respect, the noise appearing when one simulates with this method the mea- 
surement of a given observable A is also interesting. The fluctuations in the 
number of occurrences of a given eigenvalue a, of A correspond to the quantum 
noise that one would get in a real experiment, performing the relevant detection 
scheme on an individual quantum system. Since more and more quantum optics 
and atomic physics experiments are now performed with a single system (sin- 
gle ion or atom, single mode of a cavity), Monte-Carlo wave function methods 
should therefore have many applications, since they lead to predictions closer to 
actual experimental signals than the master equation, which rather deals with 
ensemble averages. 
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